Exact Numerical Study of Pair Formation with Imbalanced Fermion Populations 
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We present an exact Quantum Monte Carlo study of the attractive 1-dimensional Hubbard model 
with imbalanced fermion population. The pair-pair correlation function, which decays monoton- 
ically in the absence of polarization P, develops oscillations when P is nonzero, characteristic of 
Fulde-Ferrell-Larkin-Ovchinnikov phase. The pair momentum distribution peaks at a momentum 
equal to the difference in the Fermi momenta. At strong coupling, the minority and majority mo- 
mentum distributions are shown to be deformed, reflecting the presence of the other species, and its 
Fermi surface. The FFLO oscillations survive the presence of a confining potential, and the local 
polarization at the trap center exhibits a marked dip, similar to that observed experimentally. 
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Following the observation of transitions into Mott in- 
sulating phases as the ratio of interaction strength to 
kinetic energy is varied^, ultracold atoms trapped on 
optical lattices have been increasingly used to emulate 
the rich physics of strongly correlated condensed matter 
systems. One of the subtle phenomena being sought is 
the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase^ii^, 
in which an imbalance in the populations of up and down 
spin electrons in a superconductor leads to Bose-Einstcin 
condensation in states at non-zero momentum and to 
spatial inhomogeneities of the pairing and spin corre- 
lations. The observation of the FFLO phase in solids 
proved very difficult and was only achieved recently in 
heavy fermion systems^. 

Cold atom experiments, in which two hypcrfine states 
of fermionic atoms play the role of up and down spins, 
have now reported the presence of pairing in the case of 
unequal populations^!^. Many theoretical studies using 
mean HeldJ-^^'^^^^^^^^dM^^l^^i^, effective Lagrangianil 
and Bethe ansatai^ have been performed for the uni- 
form system with extensions to the trapped system using 
the local density approximation (LDA). For the uniform 
case, the debate revolves on the details for the paired 
state: FFLO pairs forming with non-zero momentum vs 
breached pairing (BP) at zero momentu m^^i^^i^^ , and on 
the fragility of such phases. No general consensus has 
emerged. 

In this paper wc report an exact Quantum Monte Carlo 
(QMC) study of pairing in one dimensional fermion sys- 
tems with population imbalance. We focus first on the 
homogeneous (untrapped) case since even here there is 
no agreement on the pairing mode. Our key result is 
that the FFLO phase is very robust and appears to be 
the dominant pairing mechanism. We show that the pair 
Green function exhibits oscillations characteristic of this 
phase, leading the pair momentum distribution function 
to peak at a momentum corresponding to the difference 
of the Fermi momenta of the individual species. 

In the last section, wc consider the trapped system 
and how the basic FFLO phenomena survive the result- 



ing density and polarization inhomogeneity. In addition, 
we show that at large \U\ the difference in density be- 
tween the two species exhibits a deep minimum in the 
trap center: Not only is the minority species attracted to 
the trap center by the majority potential, but in fact the 
two populations attempt to equalize there. This is a key 
signature of the experiments on unequal populations^. 

We first describe the model and the QMC method we 
employed. The attractive Hubbard Hamiltonian is, 

la 

+ C/^nnni2 + V^TX]^^("'i+'^'2), (1) 
I I 

where Cj^ricja) are fermion creation (destruction) opera- 
tors on spatial site j with the fermionic species labeled by 
a = 1,2 and nj „ = Cj ^cj a the corresponding number op- 
erator. We take the hopping parameter t = 1 to set the 
energy scale and attractive on-site interactions U < Q. 
Vt is the strength of the quadratic confining potential. 
All our results are for inverse temperature (3 = 64, so that 
T = W/256 with W = 4t the bandwidth. We have veri- 
fied that this yields ground state properties. We study a 
one dimensional lattice with L = 32 sites, unless other- 
wise stated. 

For our simulations, wc use a continuous imaginary 
time canonical "worm" algorithm where the total num- 
ber of particles is maintained strictly constan1>2i. In this 
algorithm, two worms are propagated, one for each type 
of fermion, which allows the calculation of the real-space 
Green fimctions of the two fermionic species, Ga, and 
also the pair Green function, Gpair, which are defined by, 

Gail) = 
Gpair(0 = (A,+,Aj), 

Aj = Cj2Cjl. (2) 

It is then immediate to obtain their Fourier transforms, 
riaik) and npair(fc). We emphasize that this algorithm is 
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FIG. 1: (color online). The pair Green function Gpairdi — j\) 
for iVi = 15 and iV2 = 15, 19, 23, 27 (Polarizations P = 
0,0.118,0.211, and 0.286) and U = -8. For zero polariza- 
tion, Gpair(|i — j\) decays monotonically while for P 7^ 0, it 
oscillates, indicating the presence of an FFLO state. 



exact: There are no approximations and the only errors 
are statistical (which are in all cases smaller than the 
symbol size). The polarization is given by P = {N2 — 
Ni)/ {N2 + Ni), where Ni{N2) is the minority (majority) 
particle numbers. The associated Fermi wavevectors are 
kFa = 2ttNJL. 

We begin with the uniform system, Vp = 0. In one di- 
mension, spin, charge, and pair correlations of the Hub- 
bard Hamiltonian in the ground state decay algebraically 
with increasing separation, with smallest exponent char- 
acterizing the phases^S. Here, with only an attractive 
on-site interaction, the dominant order is in the s-wave 
superconducting channel. In Fig. [T] we show the real- 
space pair correlation function Gpair(0 fo'" fixed U = —8 



and different P. Gpair(0 decays monotonically for P — 0, 
but develops clear oscillations for P 7^ 0. These oscilla- 
tions are characteristic of the FFLO state: The mismatch 
of the two Fermi surfaces leads to pairing at nonzero cen- 
ter of mass momentum and, consequently, spatially inho- 
mogcnous regions in which the pairing amplitude oscil- 
lates. The larger the Fermi surface mismatch, i.e. the 
larger the P, the larger the center of mass momentum 
and the smaller the period, as seen in Fig. [TJ The modu- 
lations in Gpair(0 follow the Larkin-Ovchinnoko^*^ (LO) 
form where the order parameter is modulated with a 
cos(gr) as a function of position r with q = ±\kpi — kp2\. 
The period, T, of the oscillations is given by T = STr/jgl 
as can be seen easily in Fig. [TJ 

To put the FFLO behavior in better context, we begin 
in Fig.[2]with n„{k) and r7,pair(fc) for the unpolarizcd case. 
At weak coupling, U = — 2, the fcrmion momentum dis- 
tribution function, 7ii(fc) ~ 71.2 (fc) has a sharp Fermi sur- 
face which is then increasingly rounded as |C/| increases. 
In all cases, n„ is a monotonic function of k. Meanwhile 
the pair momentum distribution function, npair(fc), has a 
peak at fc which grows with \U\: The pairs have zero 
momentum and become more tightly bound and numer- 



c«Dn^(k), U=-2 
(k), U=-2 

pair^ " 

oGn^(k), U=-4 
i^^n (k), U=-4 

pair^ " 

G-on^(k), U=-8 
a-aH .(k), U=-8 

pair^ " 




FIG. 2: (color online), ni(fc) and npaii(fc) for the case of 
^ N2 = 15 and L = 32, /3 = 64, (7 = -2, -4, -8. Here 
P = 0, so ni(fc) = n2{k). 
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FIG. 3: (color online), ni(fc), n2{k) and npair(fc) for the case 
of Ni = 7, N2 = 9, L ^ 32, /3 = 64 and [/ = -4 (left panel) 
and U — —10 (right panel). 



ous the larger the on-site attraction. 

When the system is polarized, it is thought to de- 
velop either pairs with non-zero momentum via the 
FFLO mechanism, or zero momentum pairs via the BP 
mechanisroi^ii^ii^. We show in Fig. [3] the single particle 
and pair momentum distributions for iVi = 7, = 9 
(P = 0.125) with J7 = -4 and U = -10 on a system 
with L — 32 sites. We note the following features: (a) 
in both cases, npair(fc) peaks at ±|A:fi — kF2\ (we show 
only fc > since the figure is symmetric), (b) the height 
of the peak increases with increasing \U\, (c) in both 
cases, the Fermi surfaces for the minority and majority 
populations are much more sharply defined than their 
P = counterparts at the same \U\. Another strik- 
ing feature in the figure is the behavior of n2{k). For 
U = —4, n2{k) displays a dip at fc = 0.39 < kF2 which 
deepens as \U\ increases and spreads to fc = 0: For ex- 
ample, for U — —10, we see that ri2(fc < 0.59) « 0.54 but 
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FIG. 4: (color online). Kinetic energy/site for the minority 
and majority populations and for the pairs as functions of \U\ 
for L = 32 and /? = 64. Also shown is the "double occupancy" 
{ni\ni2) which gives the interaction energy when multiplied 
by [/. 




FIG. 5: (color online), npair(fc) for several values of the po- 
larization showing the FFLO peak going to higher values of 
fcpoak as I^Fi — kF2\ increases. The inset shows the position 
of the peak versus \kFi — kF2\- 



then rises to n2(0.59 < k < 0.78) « 0.64 before it drops 
at the Fermi surface. This remarkable feature, in which 
the momentum distribution can increase as k increases, 
is very robust for large negative \U\ and is not a finite 
size effect. 

Similar behavior has been seen in mean field both for 
FFLO and BP in Ref. [H where the Fulde-Ferrell (FF) 
single plane wave ansatz, A(r)e"''', was used. It was 
found that in the BP case, n2{k) is symmetric around 
fc = whereas in the FFLO case n2{k) is asymmetric. 
It is likely that the LO ansatz with cosine modulation 
due to the superposition of two plane waves with oppo- 
site wavevectors will lead to a symmetric result for n2(fc) 
in the FFLO state. However, there is no ambiguity in 
the pair momentum distribution, npair(fc): its peak lies 
at fcpeak = ±|fc_Fi " kp2\ as is seen in Fig. [3] where we 
only display fc > since the figure is symmetric. Fur- 
thermore, we note that this peak starts to form at fc ^ 
even for the smallest values of |[/| we have examined. For 
example, for the case of Fig. [31 it is present aiU = —0.5. 

Additional insight on the pairing is obtained from 
the energetics. In Fig. H] we show the kinetic energy 
(KE) per site for the minority and majority populations, 
(cj+i ct'-Zct)' ^"^^ ^'^'^ pairs, (A^^^ A]^). We also show the 
subtracted pair KE/site ((A;_,_-^A;) — (c;iC;;^)(c;2C;2)) 
emphasize the contribution from pairing. We see that 
as \U\ increases, the single fcrmion kinetic energies de- 
crease (in absolute value) while the pair KE increases: 
more and more of the fcrmion hopping is performed in 
pairs. We also see clear cross-over behavior in the KE 
curves (change in the sign of the curvature) as \U\ is in- 
creased which can be understood by examining the dou- 
ble occupancy, {niini2). Clearly, (niina) lies between 
NiN2/LP' at f/ = (no pairing) and Ni/L at very large 
negative U (maximum pairing). When pairing is satu- 
rated, the system is made of A^i tightly bound pairs and 



a gas of — Ni unpaired fcrmions. The crossover in 
the KE appears to take place as the number of tightly 
bound pairs most rapidly approaches its saturation value 
a±U ~ —2.5. We emphasize that for all values of C/, the 
peak in ripair is at fcpcak 7^ when P ^ 0. 

To underscore the robustness of FFLO pairing, we 
show npair(fc) in Fig. \5\ at U — —9 for several polar- 
izations. The curves clearly show a peak at fcpcak 
which, as is seen in the inset, corresponds to {kpi — fcF2|- 
It is striking that even at the largest polarizations con- 
sidered, corresponding to A^i = 15 and N2 ^ 31, FFLO 
pairing is still present as seen in the npair(fc) peak. 
For these parameters, we have not found the Clogston- 
Chandrasekhar limii*2E in which extreme polarization 
completely destroys the superconducting state. It is pos- 
sible that going to much lower A'^i can reach it. 

Several measurements have been reported on spin im- 
balanced cold atom systems^i^ confined in 3-dimensional, 
but highly elongated, traps. In these experiments, the 
effect of the confining potential is critical to the results. 
We now apply a trap, Vt ^ 0, in our simulations and dis- 
cuss its effect on the FFLO pairing state. It is clear that 
the presence of the majority species in the trap center, 
combined with the attractive interaction, will lead to in- 
creased localization of the minority species in the center. 
Experiments^!^ show additional interesting features: this 
localization tendency is so marked that the local density 
difference ni\ — 71^2 vanishes in an extended region about 
the trap center. One recent DMRG calculation^^, ob- 
serves both the FFLO state and the squeezing of the 
minority population, but their local density difference 
is maximal at the trap center. More precisely, there 
is an extended flat region of constant local polarization 
near the center, which then falls off as the distance in- 
creases. The constancy of this polarization is reflected 
in the fact that the period of the FFLO oscillations is 
uniform throughout the central region. 
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FIG. 6; (color online). Left: Density profiles for the minority 
and majority populations in a trap Vt = 0.005. Here L — 
40, [/ = — 10, /3 = 64. The phase in the trap center is not a 
fully packed Fock state, yet the difference in the local density 
still exhibits a deep minimum. Right: The Fourier transform 
npa,ii{k) of the pair correlation function. The peak position 
is at fc = it/10, which is very close to the value k — 2n{ni — 
712) = 0.31 obtained by using the the local polarization pi — 
nn — ni2 — 0.049 at the trap center i = 20. 

A second DMRG treatment^i reports zero density dif- 
ference in the trap center associated with a Fock state 
of sites which are fully occupied with nn = ni2 = 1. 
This leads to a vanishing local polarization, but does not 
support superconducivity (at any wavevector) since the 
particles cannot move. At larger global polarization the 
Fock state melts and the minimum in local polarization 
is nonvanishing at the trap center, but the state there is 
metallic rather than superconducting, since the sites are 
still fully occupied in the majority channel. 



In Fig. [6] we show that a marked minimum in the den- 
sity difference is present even when the trap center is 
superfluid, with both nn < 1 and 71^2 < 1, as in the 
experimental situation. The figure also shows a peak in 
'T-pair(fc) at nonzero k, demonstrating that pairing occurs 
and, furthermore, that this superfluid is of the FFLO va- 
riety, as expected from the nonzero local polarization, pi. 
At the trap center, pi = nn — ni2 = 0.049 which would 
lead to FFLO pairing at ka — 0.31 in the uniform case. 
Towards the edge of the central polarization minimum, 
the maximal pi — 0.115, with an associated fcf, = 0.72. 
The observed peak in rtpair(^) lies very close to ka and is 
determined by the lower polarization region. 

In conclusion, we have studied the one-dimensional at- 
tractive Hubbard model with imbalanced populations for 
a range of P and U values with and without a trapping 
potential. In the absence of a trap, we have presented 
results for pair Green functions, single particle and pair 
momentum distributions and kinetic and interaction en- 
ergies. Our results show that FFLO pairing is surpris- 
ingly robust from very small values of |J7| all the way up 
to saturation at very large \U\ and for a wide range of P. 
We have found no values for U and P for which BP dom- 
inates over FFLO: the maximum of npair(fc) is always at 
fcpcak = ij^Fi — kp2\- We have also shown that at large 
C/|, n„{k) of the majority populations is not monotonic. 
Finally wc have shown that these features arc robust to 
the spatial inhomogeneities induced by the presence of a 
trap. 
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